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pLi Abstract 

Dynamical energy analysis was recently introduced as a new method for determining the dis- 
tribution of mechanical and acoustic wave energy in complex built up structures. The technique 
interpolates between standard statistical energy analysis and full ray tracing, containing both of 
these methods as limiting cases. As such the applicability of the method is wide ranging and addi- 
tionally includes the numerical modelling of problems in optics and more generally of linear wave 
problems in electromagnetics. In this work we consider a new approach to the method with en- 
hanced versatility, enabling three-dimensional problems to be handled in a straightforward manner. 
The main challenge is the high dimensionality of the problem: we determine the wave energy density 
both as a function of the spatial coordinate and momentum (or direction) space. The momentum 
variables are expressed in separable (polar) coordinates facilitating the use of products of univariate 
basis expansions. However this is not the case for the spatial argument and so we propose to make 
use of automated mesh generating routines to both localise the approximation, allowing quadrature 
costs to be kept moderate, and give versatility in the code for different geometric configurations. 



^ 1 Introduction 



Predicting the wave energy distribution of the vibro-acoustic response of a complex mechanical system 
is a challenging task, especially in the mid-to-high frequency regime. Standard numerical tools such as 
finite element methods become inefficient, and ray or thermodynamic approaches are often employed 
5^ . to model the wave energy fiow through the structure. Popular methods are Statistical Energy Analysis 
(SEA) [H O [3], in which the mean energy flow between subsystems is assumed to be proportional to the 
energy gradient, and the ray tracing technique, in which the wave intensity distribution is determined by 
summing over contributions of a potentially large number of ray paths |H El E]. 

SEA is in fact a low resolution ray tracing method [3 E] leading to small numerical models compared 
to ray tracing. This efficiency saving comes at a price, however: SEA has no spatial resolution of 
the energy distribution within subsystems and becomes unreliable whenever long range correlations 
in the ray dynamics are present. The recently developed Dynamical Energy Analysis (DEA) [H [9] 
provides a tool which interpolates between SEA and a full ray tracing analysis and can overcome some 
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of the problems mentioned above at a relatively small computational overhead. DEA thus enhances the 
range of applicability of standard SEA and gives bounds on the range of applicability of SEA. Related 
methods have been discussed previously in the context of wave chaos [10] and structural dynamics [TT] . 
In particular Langley's Wave Intensity Analysis (WIA) [121 113] and Le Bot's thermodynamical high 
frequency boundary element method [HI [151 [IB] include details of the underlying ray dynamics. The 
approach employed here differs from these methods by considering multiple reflections in terms of linear 
operators. Representing these operators in terms of basis function expansions then leads to SEA-type 
equations. 

In this work we develop a new approach to DEA suitable for modelling three-dimensional problems. 
The present DEA methods rely on the fact that one can easily parametrise the boundary of the region 
being modelled, and then apply an orthonormal basis approximation over the resulting boundary phase 
space coordinate system. In two dimensions this is simple as the boundary may be parametrised along 
its arc-length and the associated momentum (or direction) coordinate taken tangential to the boundary. 
The basis can be any suitable (scaled) univariate basis in both position and momentum, such as a Fourier 
basis [8] or Chebyshev polynomials [9]. Defining a suitable parametrisation for the spatial coordinate in 
three-dimensions becomes much more difficult. In momentum space spherical polar coordinates may be 
employed and so these problems do not arise. 

In order to develop a flexible code we employ automated mesh generating routines to provide a 
widely applicable parametrisation of the boundary surface for general three-dimensional structures via 
triangulation. The precision of the spatial approximation may then be improved by reflning the mesh, 
avoiding the issue of flnding a suitable basis. One avenue for potential future study stems from the 
fact that it is possible to deflne an orthogonal basis on a general triangle which reduces to Legendre 
polynomials along one edge of the domain triangle [T7]. However, in this work we restrict to a piecewise 
constant approximation on each element of the mesh for reasons of both simplicity and to keep the 
associated quadrature costs moderate for the three dimensional case. 

For the choice of momentum basis we may take a product univariate basis as mentioned above. It is 
preferable if this basis is orthogonal with respect to the standard inner product for consistency with 
both the piecewise constant spatial approximation, and the SEA limit when the lowest order momentum 
basis is applied and continuity is enforced across the mesh. The main choices are either a Fourier basis or 
Legendre polynomials. In this work we choose Legendre polynomials due to better convergence properties 
in the absence of periodic boundary conditions [18] and for consistency with the approach in [IT] should 
we wish to include a spatial basis in future work. 

The remainder of the paper is structured as follows. In Section [2l the ray tracing approximation 
is discussed and related to the Green function using short wavelength asymptotics. In Section [31 the 
concept of phase-space operators is introduced in order to represent the propagation of ray densities in 
terms of boundary integrals. The discretization of the method using spatial meshing procedures and basis 
function approximations in direction space is then detailed. Decomposition of the method for problems 
with multiple subsystems is then discussed along with links between the method and SEA. In Section [H 
the application of boundary element DEA to two-dimensional examples is discussed and verifled against 
previous work. Finally some three-dimensional examples are considered. 

2 Wave equations and asymptotics 

It is assumed that the system as a whole is characterized by a linear wave equation describing the overall 
wave dynamics including damping and radiation in a flnite domain Q G M.'^, d = 2 or 3. In this work 
only stationary problems with continuous, monochromatic energy sources are considered. We split the 
system into Nq subsystems and consider the scalar wave equation for acoustic pressure waves in each 
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homogeneous sub-domain Qi, with local wave velocity q, i = 1, Nq and Q = IJi^i Extensions to 
more complicated systems with different wave operators in different parts of the system can be treated 
with the same techniques as long as the underlying wave equations are linear, see the discussion in Ref. 

The general problem of determining the response of a system to external forcing with angular fre- 
quency w at a source point vq G can then be reduced to solving 

{k^ - H)G{r, ro; u) = -F, 5{r - tq), i = l, Nn, (1) 

with H = —A. The Green function G represents an acoustic pressure wave and Fq is a forcing amplitude 
term with units kg which is just set to unity for simplicity. The solution point is denoted r G fij and 
6 is the Dirac delta distribution. Furthermore, fcj = cu/cj + z/ij/2 is a complex valued wavenumber, where 
the imaginary part represents a subsystem dependent damping coefficient /Xj. Throughout this work we 
take i = unless used as a subscript, in which case it is an index over the number of subsystems. 

The wave energy density induced by the source is then given as 

|G(r,ro;w)P , . 

e{r, ro; u) = ^ , (2) 

for r G f2j where Qi is the density of the medium in fij. The linear wave operator H can naturally be 
associated with the underlying ray dynamics via the Eikonal approximation; for more detailed derivation, 
see Ref. [SI HSl ISO]- Using small wavelength asymptotics, the Green function in equation ([T]) may be 
written as a sum over all classical rays from ro to r for fixed kinetic energy of the hypothetical ray 
particle. One obtains |20l [21] 

Gir, ro; c) = j^-^ E A,e^''^'^~-^-^'\ (3) 



3-ro- 



where Lj is the length of the ray trajectory between ro and r including possible reflections on boundaries. 
The amplitudes Aj may be written as a product of three terms as in Ref. [S] due to damping, mode 
conversion and reflection/transmission coefficients, and geometrical factors. The phase index Uj contains 
contributions from the reflection/transmission coefficients at interfaces between subsystems and from 
caustics along the ray path. 

Analogous representations to ([3]) have been considered in detail in quantum mechanics plj and are 
also valid for general wave equations in elasticity, see Ref. for an overview. In the latter case 

G becomes matrix valued. Note that the summation in equation ([3]) is typically over infinitely many 
terms, where the number of contributing rays increases (in general) exponentially with the length of the 
trajectories included. This gives rise to convergence issues, especially in the case of low or no damping 

m- 

The wave energy density can now be expressed as a double sum over classical trajectories and 
hence 

£(r,ro;w) = C AjAj,e'''^^^^-^o'^-^^''^-'','^^l^ 

j,j':ro^r (4) 

= G [p(r, ro; u) + off-diagonal terms], 

with G = tt"^ / {gicf{27!')^'^~^^^). The dominant contributions to the double sum arise from terms in which 
the phases cancel exactly; one thus splits the calculation into a diagonal part 

p{r,ro;uj)= ^ |A,f (5) 
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where j = j', and an off-diagonal part. The diagonal contribution gives a smooth background signal 
and the off-diagonal terms give rise to ffuctuations on the scale of the wavelength. The phases related 
to different trajectories are (largely) uncorrelated and the resulting net contributions to the off-diagonal 
part are in general small compared to the smooth part, especially when averaging over frequency intervals 
of a few wavenumbers. 

It has been shown in Ref. [8j that calculating the smooth diagonal part ([5]) is equivalent to a ray 
tracing treatment. That is, the smooth part of the energy density can be described in terms of the flow 
of fictitious non-interacting particles emerging from the source point ro uniformly in all directions and 
propagating along ray trajectories. This makes it possible to relate wave energy transport with classical 
flow equations and thus thermodynamical concepts, which are at the heart of an SEA treatment. In DEA 
the classical flow is expressed in terms of linear phase space operators as detailed in the next section. 



3 Boundary Integral Formulation 

3.1 Phase Space Boundary Integral Formulation 

Following a purely kinetic viewpoint based on the interpretation that rays are trajectories of particles 
following Hamiltonian dynamics as detailed in Section 2 of Ref. [19], the time dependence of a density 
of ray trajectories (or particles) p is known to satisfy the Liouville equation 

^{X,T) + ^.Vx{p{X,r)) = 0, (6) 

where X = (r, p) denotes the phase space coordinate with position r and momentum p. The propagator 
for the Liouville equation is given by K'^{X, Y) = 6{X — Lp'^{Y)) and is the kernel of a linear phase space 
operator known as a Perron- Frobenius operator in dynamical systems theory pOt [22]. The phase space 
flow ^p'^{Y) gives the position of the particle after time r starting at Y = {r',p') when r = 0. Hence we 
may write 

p(X,r)= / K-{X,Y)po{Y)dY, (7) 



where po denotes the initial ray density at time r = 0. The domain of integration is over the whole 
of phase space P = i7 x R'^, where the integration over R"' takes care of the momentum coordinates p. 
Note that the flow satisfies the Hamilton equations of motion given by the system of ordinary differential 
equations (ODEs) 

dX_/' 1 
dr ~ [ -1 



where H = |pp is the Hamilton function for the wave operator if in ([T]). That is to say, substituting 
ip^{Y) for X(r) in (E]) satisfies the system of ODEs with X(0) = Y. 

Consider a source localized at a point ro emitting waves continuously at a fixed angular frequency 
u. Standard ray tracing techniques estimate the wave energy at a receiver point r by determining the 
density of rays starting at ro and reaching r after some unspecified time. This may be written in the 
form 

POO 1^ r 

p{r,ro,oj)= / wiY,T)K^X,Y)po{Y,u;)dYdpdT, (9) 

Jo JR'' J¥ 

with initial density pq{Y,uj) = 6{r' — ro)5(fco — H{Y)), where /cq is the wave number at the source point 
as defined in Eqn. ([T]). It can be shown that equation is equivalent to the diagonal approximation 
([5]) [8]. A weight function w is included to incorporate damping and reflection/transmission coefficients. 
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It is assumed that w is multiplicative, that is, {w{ip'^^{X),T2)w{X,Ti) = w{X,ti + r2)), which holds for 
standard absorbtion mechanism and reflection processes |20] . 

In order to solve the stationary flow problem we may rewrite equation ([9]) in boundary integral form 
using a boundary mapping technique. For the time being let us consider a problem with a single (sub- 
)system f2 = f2i with boundary F. The boundary mapping procedure involves first mapping the ray 
density emanating continuously from the source onto the boundary F. The resulting boundary layer 
density p[?'* is equivalent to a source density on the boundary producing the same ray field in the interior 
as the original source field after one reflection. Secondly, densities on the boundary are mapped back 
onto the boundary by a boundary operator B with kernel Kr{X^, Y^; u) = w{Y'^)6{X^ — (j)^{Y^)), where 
X^ = {s,Ps) represents the coordinates on the boundary. That is, s parameterizes F and Ps € B^^^ 

denotes the momentum component tangential to F at s for fixed H{X) = where B'^^^ is an open 
ball in M.'^^^ of radius \p\ and centre s. Likewise, Y'^ = {s',p'g) and 0"^ is the invertible boundary map. 
Note that convexity is assumed to ensure (p^ is well defined; non-convex regions can be handled by 
introducing a cut-off function in the shadow zone as in Ref. [15] or by subdividing the regions further. 

The stationary density on the boundary induced by the initial boundary distribution p^{X'^,uj) can 
then be obtained using 

oo 

prM = ^"(^)PrM = (/ - Biuj))-'pUoo), (10) 

n=0 

where ^B" contains trajectories undergoing n reflections at the boundary. The resulting density distri- 
bution on the boundary p^lX^^ju) can then be mapped back into the interior region. One obtains the 
density after projecting down onto coordinate space. 



3.2 Discretization and Basis Representation 

The long term dynamics are thus contained in the operator (/ — B)~^ and standard properties of Perron- 
Frobenius operators ensure that the sum over n in equation ([TU]) converges for non-vanishing dissipation, 
or in open systems. In order to evaluate (/ — B)~^ a finite dimensional approximation of the operator B 
must be constructed. In Ref. [H |9] basis expansions have been applied in both position and momentum 
coordinates, which is straightforward to implement using univariate expansions in each argument for 
i7 C M^. However, it is not straightforward to construct a general orthogonal basis with independent 
spatial arguments when f2 C M^. For this reason we employ a boundary element triangulation of F, 
with a zero order basis approximation on each element for any L^— orthonormal basis, which essentially 
results in a scaled piecewise constant boundary element approximation. This type of approximation is 
also often referred to as Ulam's method [22], although here such an approximation would be performed 
in full phase space, rather than just in its spatial component. 

For the approximation in the momentum argument we choose a basis orthogonal in for consistency 
with the spatial approximation. We choose a Legendre polynomial basis for this purpose due to good 
convergence properties without requiring periodic boundary conditions [H]. Note that for Q C M^, then 
Ps G (— and for C then ps G [0, \ki\) x [—7i,n). Denote by Ps a re-scaling of Ps to (—1, 1) 
or [—1, 1) X [—1, 1) for the two and three-dimensional cases, respectively. Let us also denote 

PUPs) = -^PmiPs) (11) 

for i7 C M^, where Pm is the Legendre polynomial of order m. For i7 C M^, m = {1711,1712) is a multi- index 
of non-negative integers. Let us write ps = {pk,Pe) and likewise Ps = {Pk,Pe)- Denote 

PUPs) = -^PmAPk)Pm,{Pe). (12) 
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Explicitly the overall approximation is then of the form 



TV 



p^(X^c^)^5^5^a^6^(^5)^m(Ps), (13) 



1=1 m=0 



where N is the order of the basis expansion, n is the number of elements and bi denotes the scaled 
(for orthonormality in an inner product) piecewise constant boundary element basis function bi{s) = 
2'^^''^^!'^ I \fAi for s in element /, and zero elsewhere. Here Ai is the surface area of element / in the 
three-dimensional case and the length of element I in the two dimensional case for I = l,..,n. Note 
that in the three-dimensional case the sum over m is a double sum over the multi-index. An analogous 
approximation is also made for pp and the values of pi^m are to be determined by solving the resulting 
linear system. 

The matrix approximation B of B{uj) for the case C is therefore given by 

Bm+l+N{l-l),l3+l+N{a~-l) = 

(2m + r 



(2m + 1) 



P^{ps)hi{s)Kr{X\Y'-uj)Pp{p'^)h^{s')dY'dX' 

w{Y^)PMp{Y')M<l>:{Yn)P^{p:)ba{s')dY'. 

Here we write (p'^ = (0^, 0^), to denote the splitting of the position and momentum parts of the boundary 
map. Also d¥ = F x is the phase space on the boundary at fixed "energy" H{X) = The only 
changes for the three-dimensional case are that the indexing is slightly more complicated due to m and 
f3 becoming multi-indices (hence the pre-factor becomes (2mi + l)(2m2 + 1)/1Q) and the definition of 
Pm changes from (fTTI) to (fT2l) . Obtaining the boundary map is not always straightforward, particularly 
for general three-dimensional geometries, and hence we write the operator in terms of trajectories with 
fixed start and end points, s' and s, as follows 



Bm+1+N{l~l),l3+1+N{a~l) 



(2m + 1) 




wiY^)PmiPsis, s'))biis)PMs, S'))bais') 



r Jr 



dp's 



ds 



(15) 

ds'ds. 



The resulting boundary integral formulation containing a pair of integrals over boundary coordinates 
bears a resemblance to standard variational Galerkin boundary integral formulations such as in 



3.3 Subsystems and links to SEA 

Recall the splitting into subsystems fij, i = 1, .., Nq introduced earlier. The dynamics in each subsystem 
are considered separately so that both variability in the wave velocity q and non-convex domains may 
be handled easily. Coupling between sub-elements can then be treated as losses in one subsystem and 
source terms in another. Typical subsystem interfaces are surfaces of reflection/ transmission due to 
sudden changes in material parameters or local boundary conditions. We describe the full dynamics in 
terms of subsystem boundary operators Bij; flow between and fli is only possible if fli n 7^ and 
one obtains an integral kernel for Bij of the form 

K,,{X:,X]) = w,,{Xt)5{Xt - C,(^I)), (16) 

where 0^- is the boundary map in subsystem j mapped onto the boundary of the adjacent subsystem i and 
X^ are the boundary coordinates of VLi. The weight Wij contains reflection and transmission coefficients 
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Figure 1: Coupled two-domain systems: configurations A, B and C respectively. 



characterizing the coupling at the interface between f2j and It also contains a damping factor of the 
form exp(— /ijL) where Hi is the damping coefficient in fij as before and L is the length of the trajectory 
from s' to s. 

Repeating the steps in the previous subsection but instead using the operator above results in a basis 
function representation, see [H |9] for details of this process. Here we employ a boundary mesh ensuring 
that the boundary of an interface between two subsystems only intersects element boundaries and not 
their interiors. An SEA treatment emerges when approximating the individual operators Bij in terms of 
constant functions only [8]. Here this corresponds to an approximation in terms of the lowest order basis 
functions only, together with a coarse spatial mesh with only one element per subsystem, or more typically 
a piecewise constant approximation on a mesh with continuity enforced within each subsystem. In this 
case the matrix B can be reduced to one element per subsystem interaction, say between subsystem i and 
subsystem j, possibly with i = j. The matrix entry gives the mean transmission rate from subsystem 
j to subsystem i. It is thus equivalent to the coupling loss factor used in standard SEA equations [2]. 
The resulting full A'^q- dimensional B matrix yields a set of SEA equations using the relation f[T(I|) after 
mapping the boundary densities back into the interior. 

4 Numerical results 

4.1 Verification in 2D for Coupled Two- Cavity Problems 

In this section we consider two-dimensional polygonal domains whose boundaries are meshed by subdi- 
viding each side into equidistant sections governed by a mesh parameter Ax. The number of elements 
on any given side is computed using the integer part of the side length divided by Ax. The Jacobian 
from ( 1T5|) is written in the form 



where n and n' are the internal unit normal vectors to F at r and r', respectively. In order to treat 
the corner singularities in ( !T7|) . Gaussian quadrature is employed where end-points are not included as 
quadrature points. The convergence of the quadrature rules is still slow due to the peak in the integrand 
at corners. Telles' transformation techniques are employed to speed up the convergence p4| . 

A number of two-cavity systems are considered as shown in Fig. [T] with Dirichlet boundary conditions 
on the outer boundary. Configuration A features irregular shaped, well separated pentagonal subsystems. 
In configuration B the size of the interface between the subsystems is increased reducing their dynamical 
separation. Configuration C includes a rectangular left-hand subsystem channelling rays out of the 
subsystem and introducing long-range correlations in the dynamics. In addition, the source is further 
from the intersection of the two subsystems. Note that SEA results are in general insensitive to the 
position of the source, whereas actual trajectory calculations may well depend on the exact position. 




(17) 



7 



20 




Frequency (Hz) 



Figure 2: (Color online) Ratio of total energies R = ||Gip/||G2p in configurations A, B and C respec- 
tively. The dotted lines correspond to FEM calculations of the full wave problem. 

In [SI [9] it is demonstrated, as expected, that SEA works well for configuration A, but not so well 
for configurations B and C. In this communication we seek to verify our new approach against results 
from previous work. In particular we discuss the relative computational efficiency of the new and old 
approaches and how they scale as the level of precision in the model is increased. Energy distributions 
have been studied as a function of the frequency with a hysteretic damping factor t] = 0.01, where 
/ij = ujri/{2ci) for i = 1,2. Here and in the remainder of this section the subsystems are numbered 1,2 
from left to right. The other parameters are set to unity for simplicity, that is Qi = Ci = 1 for i = 1,2. 
For this reason the subsystem interface reflection and transmission coefficients appearing in the weight 
term in ( IT6l) are simply and 1, respectively. 

Fig. |2] shows the the boundary element DEA results for configurations A, B and C compared with 
discontinuous Galerkin finite element method simulations as detailed in |9]. Explicitly we compute the 
energy ratios between the two subsystems R = HGi ||^/||G'2||^ where 

||G,f := / \G{r,ro;uj)\^dr, t = 1,2. (18) 

The dotted lines each represent a simulation at a different frequency in the range ±5Hz of the frequencies 
used for the boundary element DEA calculations. In all three cases good convergence of the method is 
demonstrated. For configuration B, shown in the central subplot, the results converge with a slightly 
lower order of approximation. This may be due to the irregular geometry and the wide opening linking 
the subsystems meaning that the energy is very evenly distributed throughout the whole domain. Hence 
lower order spatial approximations will be reasonably good. 
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Table 1: Computational times for Configuration A with / = lOHz 
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N 


Total Computational times (s): 
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Figure 3: Energy distribution along the x-axis in the cuboid example. 

Table [1] shows the total computational times for the lOHz calculation in Fig. |2] using both boundary 
element DEA and comparing with a previous approach where a Chebyshev basis is employed in full phase 
space [9]. The computations were performed using a desktop PC with a 2.83 GHz dual core processor, 
although the code was not parallelized. The total computational expense is considerably reduced using 
the current boundary element DEA approach. In addition the computational cost of boundary element 
DEA is growing more slowly as the precision of the model is increased. This will be very important for the 
three-dimensional case where the number of degrees of freedom in the model will increase more quickly. 
When we consider that the Chebyshev algorithm was already a considerable saving on the original DEA 
methods discussed in [8] one can see how far we have come. Since the parameters for configurations B 
and C are similar to those for configuration A, the computational times are roughly the same for the 
same orders of approximation. 

4.2 Applications in 3D 

In this section we consider some three-dimensional domains whose boundaries have been triangulated 
using the Tetgen freeware automated mesh generating package (http:\\tetgen. berlios.de). The 
Jacobian from f|T5l) may be computed using standard formulae for global, local and polar coordinate 
transformations along with several applications of the chain rule. As before the Jacobian introduces 
singularities in the integrals along edges and at vertices of the domain. Again Gaussian quadrature and 
a carefully chosen application of Telles' transformation techniques are employed to ensure relatively fast 
convergence of the numerical integration procedures. 
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Figure 4: The car cavity showing the source point and the open roof (black) acting as absorbing boundary. 

The first example we consider is that of a cuboid {x,y,z) G (—1, 1) x (—0.5,0.5) x (—0.5,0.5) with 
Dirichlet boundary conditions. The source point is taken as (—0.9,0.1,0.1) and the same frequency 
and damping correspondence is used as in the two-dimensional examples and g = c = 1. Figure [3] 
shows the computed energy distributions inside the cuboid along the x-axis. The method is compared 
against discontinuous Galerkin finite element method (FEM) computations, which are averaged over 17 
frequencies within ±2Hz of the (central) frequency used for the boundary element DEA computation. 
Further details of the FEM techniques employed here can be found in |9] and references therein. The 
dashed line shows an approximation with a coarse mesh and where the energy density is assumed constant 
over all possible directions of rays approaching the boundary from the interior. The computation time 
for such an approximation is typically around a minute per frequency. The dash-dot line shows a higher 
order approximation where we have refined the mesh until the solution appears reasonably converged by 
eye, in this case 344 elements were employed. Also due to the relatively low dissipation levels in these 
plots a low order approximation in momentum (quadratic) was sufficient to give reasonably converged 
results. The computation time for this plot was approximately 16 hours using a parallel machine with 
two quad-core processors. 

It has been demonstrated that in the low damping regime SEA can provide reasonably good ap- 
proximations even in regular structures [9], which explains why the coarse approximation here is still 
reasonably good. However, one notices an improvement in the match with the FEM data at both the 
peak and tail of the plot for both frequencies considered when the higher order approximation method is 
employed. There is, however, a significant computational cost associated with this increased precision. 

The second example we consider is an open car cavity as discussed in [25] and shown in Fig. SJ The 
source point is located on the base of the cavity at (0.6, 0.0, 0.4). Again we consider Dirichlet boundary 
conditions except along the roof, which is assumed to be non-refiecting along the subsection shown in 
black in Fig. |U between x = 1.0 and x = 1.8. Physically this corresponds to an opening in the cavity 
with identical media in both the interior and exterior. 

Figu. E] shows the amplitude |G| plotted in the interior of the cavity along the midplane z = 0.85. 
The jagged diagonal edge shown in Fig. [5] is merely an artifact of the plotting and interpolation of the 
amplitude at discrete points, and is actually smooth in the model as shown in Fig. HI No damping 
is incorporated in this example and energy losses only occur through the open roof, meaning that 
the plots here are now wavenumber independent. The three subplots show successively higher order 
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Figure 5: Amplitude \G\ along the plane z=0.85 for the car cavity example. 
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approximations from upper to lower and convergence in the plots is evident due to the increased similarity 
between the lower two plots. Directivity in the wave field plays a much stronger role in this example due 
to the localised dissipation at the opening of the cavity. For this reason it was necessary to employ a 
higher order momentum basis approximation than before with = 4 before the plot appears reasonably 
converged. The computational times were approximately one minute for the upper subplot, eleven hours 
for the central subplot and three days for the lower plot, again using a parallel machine with two quad-core 
processors. 

The upper plot in Figure |5] has a markedly different appearance to the other subplots showing that the 
most coarse approximation is not good in this case. One reason for this is the much stronger directivity 
of the wave field compared with the previous example. One can clearly see how in the upper plot the 
solution is more slowly varying and distributed more uniformly as you move away from the source point. 
In the lower plot one sees a noticeable dip in the amplitude close to the non-refiecting boundary. Also the 
increased intensity around the source point stretches more in the horizontal direction than the vertical 
direction in the lower plot, but is more evenly distributed in the other plots. In all three subplots the 
wave amplitude is greater in the region to the left of the opening (for x < 1) since in this region there are 
many possible ray trajectories that remain trapped in the cavity for a long time before exiting through 
the opening. This is also true for the region y < 0.8 and hence the greatest intensities are observed in the 
intersection of these two regions. In billiard dynamics these trajectories are known as near-bouncing-ball 
orbits, see for example [26]. 

5 Conclusions 

A new approach to determining the distribution of mechanical and acoustic wave energy in complex built 
up structures has been discussed. The methodology has been carefully chosen to permit application to two 
or three dimensional problems. Using boundary element meshes for three dimensional problems renders 
the method applicable to general domains and removes the need to determine an orthogonal spatial 
basis for each geometrically different example. The application of the method to some well studied 
two-dimensional examples has shown it to be relatively efficient and scale favourably as the number of 
degrees of freedom in the model is increased compared with previous DEA approaches. Examples in three 
dimensions were also considered showing both the applicability and versatility of the method, but also 
its high computational cost as the number of degrees of freedom is increased. We have also seen however 
that in some cases a low order and fast computation yields reasonably good results. The suitability of 
the method for parallel processing means that with greater computing resources it has the potential to 
be employed in larger and more complicated configurations than those considered here. 
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